function [dates, data] = gFilter(dates, data, minFreq, maxFreq, toPlot)
%GFILTER Apples the specified filter to the data
%   [DATES,DATA] = GFILTER(DATES,DATA,FILTERTYPE,FILTERPARAM) Filters the
%    supplied data with cutoff frequencies equal to minFreq and maxFreq.
%    
%    Cutoff frequencies should be specified in cycles per year, and dates
%    should be given in Matlab datenum format.
%
%    To create a High Pass Filter, set maxFreq to zero.
%
%    To create a Low Pass Filter, set minFreq to zero.
%
%   GFILTER(...,'plot') filters the data and, in addition, plots the
%   results as a timeseries line plot.

Ts = 365/mean(diff(dates));                        % Samples per year

xfft = fft(data);                                  % Data transformed to frequency domain 

N = length(data);

tt = [((1:ceil(N/2))-1)/N (floor(N/2):-1:1)/N];    % Vector of digital frequencies
                                                   
freqs = tt.*Ts;                                    % Create a vector of frequencies
                                                   
if(minFreq==0)                                     % Zero out frequency coefficients for frequencies outside the specified bounds
    lowerBound = [];                                
    filterType = ['Low Pass Filter with a cutoff frequency of ',num2str(maxFreq),' cycles/year'];
else                                                
    lowerBound = freqs<minFreq;
    filterType = ['Band Pass Filter with a passband of ',num2str(minFreq),' to ',num2str(maxFreq),' cycles/year'];                      
end
if(maxFreq==0)
    upperBound = [];
    if(strcmpi(filterType,['Low Pass Filter with a cutoff frequency of ',num2str(maxFreq),' cycles/year']))
        filterType = 'Unity Gain Filter';
    else
        filterType = ['High Pass Filter with a cutoff frequency of ',num2str(minFreq),' cycles/year'];
    end
else
    upperBound = freqs>maxFreq;
end
xfft(lowerBound) = 0;
xfft(upperBound) = 0;

data = ifft(xfft);                                % Transform the remaining frequency coefficients back to the time domain
                                                  
if(nargin < 5); return; end;                      

if(strcmpi(toPlot,'plot'))                        % Plot the filtered data
    figure                                        
    plot(dates,data);
    datetick('x');
    xlabel('Time');
    title(['Data after passing through a ',filterType]);
end

end

